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We  have  conducted  N-body  simulations  of  the  growth  of  Milky  Way-sized  halos  in  cold  and 
warm  dark  matter  cosmologies.  The  number  of  dark  matter  satellites  in  our  simulated  Milky  Ways 
decreases  with  decreasing  mass  of  the  dark  matter  particle.  Assuming  that  the  number  of  dark 
matter  satellites  exceeds  or  equals  the  number  of  observed  satellites  of  the  Milky  Way  we  derive 
lower  limits  on  the  dark  matter  particle  mass.  We  find  with  95%  confidence  tUs  >  11.8  keV  for  a 
sterile  neutrino  produced  by  the  Dodelson  &  Widrow  mechanism  and  mwDM  >2.1  keV  for  a  thermal 
dark  matter  particle.  The  recent  discovery  of  many  new  dark  matter  dominated  satellites  of  the 
Milky  Way  in  the  Sloan  Digital  Sky  Survey  allows  us  to  set  lower  limits  comparable  to  constraints 
from  the  complementary  methods  of  Lyman-a  forest  modeling  and  the  unresolved  cosmic  X-ray 
background.  Future  surveys  like  LSST,  DES,  PanSTARRS,  and  SkyMapper  have  the  potential  to 
discover  many  more  satellites  and  further  improve  constraints  on  the  dark  matter  particle  mass. 


I.  INTRODUCTION 

Cold  dark  matter  (CDM)  is  extremely  successful  at  de¬ 
scribing  the  large  scale  features  of  matter  distribution  in 
the  Universe  but  has  problems  on  small  scales.  Below  the 
Mpc  scale  CDM  predicts  numbers  of  satellite  galaxies  for 
Milky  Way-sized  halos  about  an  order  of  magnitude  in 
excess  of  the  number  observed.  This  is  the  ‘missing  satel¬ 
lites’  problem  Hi-  One  proposed  solution  is  that,  due 
to  feedback  mechanisms,  dark  matter  satellites  below  a 
certain  scale  do  not  form  stars  and  are  nonluminous  dark 
halos.  Another  solution  is  the  dark  matter  may  be  ‘warm’ 
(particle  mass  ~  1  keV)  instead  of  ‘cold’  (particle  mass 
~  1  GeV).  By  reducing  the  dark  matter  particle  mass, 
streaming  motions  can  erase  density  fluctuations  on  small 
scales  and  reduce  the  number  of  satellites.  Warm  dark 
matter  (WDM)  models  have  been  studied  through  N- 
body  simulations  by  a  number  of  authors  il. 

A  complication  arises  when  running  N-body  simula¬ 
tions  of  WDM  cosmologies  due  to  numerical  artifacts 
produced  by  discrete  sampling  of  the  gravitational  poten¬ 
tial  with  a  finite  number  of  particles.  That  discreteness 
effects  are  a  problem  for  simulations  started  from  lattice¬ 
like  initial  conditions  has  long  been  known  (see  @  for  a 
review).  Particles  are  positioned  on  a  regular,  uniform 
lattice  and  then  their  positions  are  perturbed  by  a  ran¬ 
dom  realization  of  the  linear  fluctuation  field  associated 
with  the  choice  of  cosmology.  Perturbations  collapse  and 
form  filaments  with  nonphysical  clumps  separated  by  a 
distance  equal  to  the  lattice  spacing  (see  Figure[T]).  Wang 
&  White  [13  have  shown  that  discreteness  also  affects 
simulations  started  from  glass-like  conditions.  A  glass 
is  created  by  advancing  a  random  distribution  of  parti¬ 
cles  under  repulsive  gravitational  forces  until  the  system 
reaches  quasi-equilibrium  and  the  total  force  on  each  par- 
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tide  vanishes.  The  particle  positions  are  then  perturbed 
from  this  state  as  for  the  lattice.  Wang  &  White  showed 
that  even  under  these  conditions,  with  no  large-scale  or¬ 
der,  filaments  still  fragment  into  halos  with  separation 
equal  to  the  mean  particle  separation  in  the  glass.  These 
halos  are  nonphysical  numerical  artifacts.  The  ability 
of  these  halos  to  survive  disruption  as  they  accrete  from 
filaments  onto  a  Milky  Way-sized  halo  has  not  been  stud¬ 
ied,  but  they  may  contaminate  the  satellite  distribution 
in  WDM  simulations. 


FIG.  1.  “Poisson  noise”  halos  formed  along  a  filament  and 
accreting  onto  a  larger  halo  at  «  =  1  in  a  WDM  simulation 
(mwDM  =  1  keV).  These  halos  are  numerical  artifacts. 

In  the  last  few  years  16  new  dwarf  spheroidal  galax¬ 
ies  have  been  discovered  in  the  Sloan  Digital  Sky  Survey 
(SDSS)  [ni  (see  Table  3  and  references  therein).  Af¬ 
ter  correcting  for  completeness  the  estimated  number  of 
Milky  Way  satellites  is  >  60.  These  new  dwarfs  have  low 
luminosities,  low  surface  brightnesses,  and  are  dark  mat¬ 
ter  dominated.  Since  the  number  of  dark  matter  halos 
must  be  greater  than  or  equal  to  the  number  of  observed 
satellites,  the  new  data  from  the  SDSS  may  provide  in¬ 
teresting  limits  on  how  cold  the  dark  matter  is. 

Motivated  by  the  recent  increase  in  the  number  of  ob¬ 
served  Milky  Way  satellites  and  by  the  need  to  better 
understand  discreteness  effects,  we  have  performed  new 
simulations  of  the  growth  of  a  Milky  Way-like  galaxy  in 
CDM  and  WDM  cosmologies  for  a  variety  of  WDM  par¬ 
ticle  masses.  Our  goal  is  to  constrain  the  dark  matter 
particle  mass  by  comparing  the  number  of  satellite  ha- 
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los  in  the  simulated  Milky  Way  to  the  observed  number 
of  luminous  satellites  for  the  actual  Milky  Way.  Maccio 
&  Fontanot  Q  combined  N-body  simulations  with  semi- 
analytic  models  of  galaxy  formation  to  compare  the  sim¬ 
ulated  and  observed  Milky  Way  satellite  luminosity  func¬ 
tions  for  CDM  and  WDM  cosmologies.  In  this  work,  we 
do  not  make  any  assumptions  on  how  we  populate  dark 
matter  halos  with  luminous  galaxies.  We  simply  impose 
that  the  number  of  observed  satellites  is  less  than  or  equal 
to  the  number  of  dark  matter  halos  for  a  range  of  Galac- 
tocentric  radii.  This  guarantees  a  robust  lower  limit  on 
the  dark  matter  particle  mass. 


II.  SIMULATIONS 


All  our  simulations  were  conducted  with  the  N-body 
cosmological  simulation  code  GADGET-2  [l2|  assuming 
dark  matter  only  simulations  (we  do  not  include  the  effect 
of  baryons  in  our  simulations).  We  adopted  values  for 
cosmological  parameters  from  the  third  year  release  of 
the  WMAP  mission  [I^,  (0^,  Da,  h,  ag,  rig)  =  (0.238, 
0.762,  0.73,  0.751,  0.951).  For  each  simulation  set  we 
produced  a  single  realization  of  the  density  field  in  the 
same  periodic,  comoving  volume  but  varied  the  power 
spectrum  of  fluctuations  appropriate  for  CDM  and  WDM 
cosmologies.  Our  initial  conditions  were  generated  on  a 
cubic  lattice  using  the  GRAFIC2  software  package  H- 
The  power  spectra  for  CDM  and  WDM  are  given  by 

PcDMik)  oc  (1) 

PwDM(k)  =  PcdmTwdm^  (2) 


respectively,  with  the  normalization  of  Pc  dm  determined 
by  cTg.  We  used  the  transfer  function  for  CDM  adiabatic 
fluctuations  given  by  Bardeen  et  al.  (BBKS)  m-- 


TcDM{k) 


ln(l  -h  2.34q) 
2.34g 


[1  -f  3.899  +  (le.lg)^  +  (5.469)3  (6.719)^]  (3) 


where  q  =  A  potential  problem  with  the 

BBKS  transfer  function  is  that  it  underestimates  power 
on  large  scales.  In  the  Appendix  we  investigate  the  ef¬ 
fect  that  our  choice  for  the  CDM  transfer  function  may 
have  on  the  number  of  Milky  Way  satellites.  We  run  one 
of  our  simulations  adopting  the  transfer  functions  from 
Eisenstein  &  Hu  M  and  we  hnd  that  this  does  not  affect 
our  results  on  the  number  of  satellites. 

Assuming  the  WDM  to  be  a  thermal  particle,  a  particle 
that  was  in  thermal  equilibrium  with  the  other  particle 
species  at  the  time  of  its  decoupling,  we  used  the  trans¬ 
fer  function  valid  for  thermal  particles  given  by  Bode, 
Ostriker,  &  Turok 


TwDM{k)  =  [1  +  (ak/h)’']  (4) 


where  v  =  2.4,  /r  =  4.167  and 


The  parameter  gx  is  the  number  of  degrees  of  freedom 
for  the  WDM  particle,  conventionally  set  to  the  value 
for  a  light  neutrino  species:  gx  =  1-5.  The  parameter  k 
is  the  spatial  wavenumber  in  Mpc~^  and  raw  dm  is  the 
mass  of  the  WDM  particle  in  keV. 

A  candidate  for  a  thermal  WDM  particle  is  the  grav- 
itino,  the  superpartner  of  the  graviton  in  supersymmetry 
theories.  The  lightest  stable  particle  (ESP)  in  supersym¬ 
metry  theories  are  natural  dark  matter  candidates.  If 
the  scale  where  supersymmetry  is  spontaneously  broken 
is  <  10®  GeV,  as  predicted  by  theories  where  supersym¬ 
metry  breaking  is  mediated  by  gauge  interactions,  then 
the  gravitino  is  likely  to  be  the  lightest  stable  particle 
and  can  have  a  mass  reaching  into  the  keV  regime. 

In  general  the  dark  matter  particle  may  not  have  been 
in  thermal  equilibrium  when  it  decoupled.  This  is  the 
case  for  a  sterile  neutrino  (see  and  references  therein) , 

a  theoretical  particle  added  to  standard  electroweak  the¬ 
ory,  the  only  matter  it  interacts  with  (except  through 
gravity)  are  left-handed  neutrinos.  There  are  several 
mechanisms  by  which  sterile  neutrinos  can  be  produced 
.  In  the  standard  mechanism  proposed  by  Dodelson  & 
Widrow  (DW)  [l^,  sterile  neutrinos  are  produced  when 
oscillations  convert  some  of  the  more  familiar  active  neu¬ 
trinos  into  the  sterile  variety.  The  amount  produced  de¬ 
pends  on  the  sterile  neutrino  mass  and  the  mixing  angle 
but  we  will  not  consider  such  details  here  and  when  con¬ 
sidering  sterile  neutrinos  we  simply  assume  they  compose 
the  entirety  of  the  dark  matter.  The  transfer  function  for 
DW  sterile  neutrinos  with  mass  mg  is  given  by  di: 


T«(A:)  =  [l  +  (afc/h)^■^  (6) 


where  a  =  2.25,  g  =  3.08  and 


a  =  0.1959 


mg  y 

■  1  keV  / 


-0.858 


0.238 


-0.136 


0.73 


0.692 


Viel  et  al.  [2^  give  a  scaling  relationship  between  the 
mass  of  a  thermal  particle  and  the  mass  of  the  DW  ster¬ 
ile  neutrino  for  which  the  transfer  functions  are  nearly 
identical: 


mg  =  4.379  keV 


mwDM^^/^  f  Dm  \  f  h  \ 
1  keV  )  VO-2387  J 


-2/3 


(8) 

Other  sterile  neutrino  production  mechanisms  include 
that  of  Shi  &  Fuller  (SF)  [2l|  who  showed  the  DW  mech¬ 
anism  is  altered  in  the  presence  of  a  universal  lepton 
asymmetry  where  production  can  be  enhanced  by  res¬ 
onance  effects.  Sterile  neutrinos  can  also  be  produced 
from  decays  of  gauge-singlet  Higgs  bosons  at  the  elec¬ 
troweak  scale  I22II.  The  momentum  distribution  of  the 
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sterile  neutrinos  depends  on  the  production  mechanism. 
In  the  absence  of  transfer  function  calculations  we  use 
the  expressions  in  [13  for  the  free  streaming  length  and 
average  momentum  to  derive  approximate  scaling  fac¬ 
tors  for  the  SF  and  Higgs  produced  sterile  neutrinos: 
rriDw/msF  =  1-5,  mow/mHiggs  =  4.5. 

There  are  ways  other  than  warm  dark  matter  to  re¬ 
duce  small  scale  power.  Broken-scale  invariance  inflation 
models  [2^  have  a  cutoff  length  below  which  power  is 
suppressed.  Particle  theories  where  the  LSP  dark  matter 
particle  arises  from  the  decay  of  the  next  lightest  super- 
symmetric  particle  (NLSP)  can  also  suppress  small  scale 
power  if  the  NLSP  is  charged  and  coupled  to  the  photon- 
baryon  plasma  [ISl  or  if  the  NLSP  decay  imparts  a  large 
velocity  to  the  LSP  [l^.  Our  method  can  be  applied  to 
constrain  these  models  as  well.  However,  in  the  rest  of 
this  paper  we  will  not  discuss  further  the  consequences 
that  our  work  has  on  these  theories. 

In  our  simulations  we  assume  the  dark  matter  is  ther¬ 
mal  and  scale  the  results  to  the  standard  sterile  neutrino 
mass  using  Eq.  ([8]).  Our  initial  conditions  include  par¬ 
ticle  velocities  due  to  the  gravitational  potential  using 
the  Zeldovich  approximation  but  we  do  not  add  random 
thermal  velocities  appropriate  for  WDM  to  the  simula¬ 
tion  particles.  Bode,  Ostriker,  &  Turok  Q  argue  that  for 
warm  particle  masses  greater  than  1  keV  thermal  mo¬ 
tions  are  unimportant  for  halos  on  scales  of  a  kiloparsec 
and  above.  Regardless,  we  expect  thermal  motions,  if 
anything,  would  reduce  the  number  of  small  mass  halos 
and  by  not  including  thermal  motions  the  mass  limits 
derived  from  our  simulations  will  be  more  conservative. 

We  ran  simulations  for  CDM  and  WDM  cosmologies 
with  WDM  particle  masses  of  mwDM  =  1,  2,  3,  4,  and 
5  keV  (tos  =  4.4,  11.0,  18.9,  27.8,  37.4  keV).  Figured 
shows  the  power  spectra  for  these  cosmologies  along  with 
the  spectrum  for  an  11  keV  standard  sterile  neutrino  us¬ 
ing  Eq.  dS]).  We  ran  two  separate  sets  of  simulations 
both  consisting  of  a  comoving  cubic  box  90  Mpc  on  a 
side.  Set  A  consisted  of  204^  particles  giving  a  ‘coarse’ 
particle  mass  of  3.0  x  10®  Mq  and  a  force  resolution  of 
8.8  kpc  (all  our  force  resolutions  were  fixed  in  coinov- 
ing  coordinates) .  We  ran  the  HOP  halo  finding  software 
at  z  =  0  and  identified  Milky  Way-sized  halos  with 
masses  1  —  2  x  10^^  Mq.  Halos  were  examined  visually, 
one  was  chosen  that  was  at  least  several  Mpc  away  from 
clusters  and  other  large  structures  so  as  to  be  relatively 
isolated.  Its  particles  were  identified  in  the  initial  con¬ 
ditions  and  a  cubic  refinement  level,  6.2  Mpc  on  a  side, 
was  placed  on  the  region.  For  the  refinement  region  in 
our  low  resolution  simulations  we  used  11,  239, 424  (224^) 
particles  with  mass  and  force  resolutions  of  7.3  x  10^  Mq 
and  550  pc,  respectively.  For  CDM  and  WDM  parti¬ 
cle  masses  of  1,  2,  and  4  keV  we  also  conducted  higher 
resolution  simulations  with  89,  915,  392  (448^)  particles 
in  the  refinement  region  and  mass  and  force  resolutions 
of  9.2  X  10"*^  Mq  and  275  pc,  respectively.  The  simu¬ 
lated  Milky  Way  halo  had  a  neighbor  halo  with  mass 
0.23Mmw  at  a  distance  of  700  kpc  in  the  low  resolution 


FIG.  2.  Power  spectra  for  our  simulations.  The  dotted  line  is 
the  power  spectrum  for  an  11  keV  standard  sterile  neutrino 
from  Abazajian  [T^ .  The  neutrino  spectrum  is  approximately 
the  same  as  a  2  keV  thermal  particle,  validating  the  scaling 
relation  of  Viel  et  al.  [^.  The  vertical  dashed  lines  indicate 
the  lattice  cell  size  in  our  high  and  low  resolution  refinement 
levels. 

simulations.  The  real  Milky  Way  has  a  massive  neighbor 
in  M31,  the  Andromeda  galaxy  (M^nd  ~  1  —  3Mmw), 
at  a  distance  ^  700  kpc.  But  in  the  higher  resolution 
simulation  this  satellite  is  merging  with  the  Milky  Way 
at  z  =  0.  Such  a  merger  may  disrupt  the  equilibrium 
of  the  halo  and  makes  it  nonrepresentative  of  the  actual 
Milky  Way.  The  difference  between  the  high  and  the  low 
resolution  simulations  is  significant  and  complicates  the 
comparison  between  the  simulations. 

The  need  to  explore  the  scatter  between  the  satellite 
distribution  of  different  realizations  of  Milky  Way-type 
halos,  in  addition  to  the  complications  arising  with  the 
high  and  low  resolution  simulations  of  set  A,  prompted 
us  to  conduct  a  new  set  of  simulations.  Set  B  con¬ 
sisted  of  408^  particles  giving  a  coarse  particle  mass  of 
3.8  X  10®  Mq  and  a  force  resolution  of  4.4  kpc.  We  ran 
HOP  and  identified  halos  with  masses  0.8  — 2.2  x  10^^  Mq. 
For  each  halo  we  also  found  the  nearest  neighboring  halo 
with  mass  >  0.8  x  10^®  Mq.  We  selected  a  halo  whose 
nearest  massive  neighbor  was  at  least  5  Mpc  away  and 
visually  verihed  the  halo  was  indeed  isolated.  A  rect¬ 
angular  refinement  level  6.1  x  7.0  x  7.9  Mpc  was  placed 
over  this  halo’s  particles  in  the  initial  conditions.  Low 
and  high  resolutions  were  conducted  with  the  same  mass 
and  force  resolutions  as  set  A.  The  low  resolution  sim¬ 
ulations  used  16,515,072  255®)  particles  in  the  re¬ 

finement  level  while  high  resolution  used  132, 120, 576 
(~  510®)  particles  in  the  refinement  level. 

Table  U  summarizes  the  properties  of  our  simulated 
Milky  Way  halos  at  2;  =  0  in  all  the  simulations.  We 
write  i?A  to  mean  the  radius  enclosing  an  overdensity  A 
times  the  critical  value.  The  mass  and  number  of  parti¬ 
cles  inside  Ra  are  Ma  and  Na,  respectively;  va  is  the 
circular  velocity  v'^  =  GMa  / Ra  at  Ra  ,  and  Vmax  is  the 
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maximum  circular  velocity  of  the  halo. 

Figure  |3]  shows  the  density  prohles  of  the  A  and  B 
Milky  Way  halos  calculated  by  breaking  the  halo  into 
spherical  shells.  Small  differences  between  the  high  and 
low  resolution  set  A  halos  caused  by  the  merging  neighbor 
are  apparent  but  generally  the  profiles  are  very  similar 
across  all  simulations  of  each  set.  We  do  not  see  an  inner 
flattening  of  the  halos  in  the  WDM  simulations  because 
we  did  not  add  thermal  motions  to  our  particles.  Fig¬ 
ure  |4]  shows  portraits  of  the  Milky  Way  in  our  set  B  low 
resolution  simulations. 


FIG.  3.  Density  profile  of  Milky  Way  halo  in  all  simulations. 
Thick  lines  are  the  high  resolution  simulations.  The  set  A 
simulations  and  the  WDM  cosmologies  in  each  set  have  been 
offset  downward  for  clarity.  The  profiles  are  plotted  start¬ 
ing  from  the  convergence  radius  of  Power  et  al.  for  both 
resolutions  (vertical  lines). 


A.  Identification  of  Satellites 

Finding  gravitationally  bound  dark  matter  halos  in  N- 
body  simulations  is  accomplished  with  halo  finding  soft¬ 
ware  employing  a  variety  of  algorithms.  Traditional  halo 
finders  like  hOF,  DENMAX  and  its  relatives  HOP  and 
SKID  [H,  [2814^.  have  difficulty  finding  satellite  halos 
within  a  much  larger  host  halo.  New  satellite-savvy  halo 
finders  such  as  6DFOF,  AHF,  and  SUBFIND  [mHII, 
have  been  developed  but  require  greater  computational 
resources  and  are  not  always  publicly  available.  Our  ap¬ 
proach  is  to  first  subtract  the  average  density  profile  of 
the  Milky  Way  from  the  particle  densities,  remove  low 
density  particles  from  the  list,  then  run  a  traditional  halo 
finder  on  the  satellite  cores  that  remain. 

We  calculated  initial  particle  densities  using  the  HOP 
software  package  smoothing  over  the  10  nearest  neigh¬ 
bors.  We  broke  the  Milky  Way  halo  into  concentric  shells 
with  equal  thickness  in  logarithmic  radius  from  0.1  to 
500  kpc.  In  order  to  account  for  the  triaxiality  of  halos, 
the  shells  were  divided  into  8  equal  volume  regions  and 
the  density  in  each  was  calculated.  We  subtracted  these 


densities  from  the  HOP  particle  densities.  The  number 
of  shells  was  adjusted  until  good  subtraction  of  the  Milky 
Way  profile  was  achieved,  see  Figured  We  also  tried  el¬ 
lipsoidal  shells  aligned  with  the  axes  of  the  halo  but  the 
improvement  was  negligible.  Particles  with  low  densities 
were  removed  and  we  ran  SKID  on  the  remaining  parti¬ 
cles  with  a  linking  length  of  2  kpc  and  smoothing  over  32 
neighbors.  Unbound  particles  were  iteratively  removed 
and  we  selected  gravitationally  bound  halos  with  10  or 
more  particles. 

SKID  calculates  properties  of  the  halos  it  finds  such  as 
the  total  mass  and  the  maximum  circular  velocity.  For 
our  purposes  the  maximum  circular  velocity  is  a  better 
characteristic  of  a  satellite  than  the  mass  because  quan¬ 
tifying  the  outer  boundary  of  a  satellite  embedded  in  a 
larger  halo  is  complicated.  Our  subtraction  of  the  Milky 
Way  profile  to  isolate  the  satellites  may  also  affect  their 
outer  boundaries.  The  max  circular  velocity  however  typ¬ 
ically  occurs  at  a  radius  well  inside  the  satellite  outskirts. 


III.  RESULTS 

A.  Satellite  Distribution  Functions 

Figure  |6]  shows  the  cumulative  max  circular  velocity 
function,  ]V(>  Vmax),  for  all  satellites  within  R^q  of  each 
simulated  Milky  Way.  Also  plotted  are  the  fitting  func¬ 
tions  of  N{>  Vmax)  from  the  Via  Lactea  I  and  H  and 
Aquarius  simulations  [H,  [H,  [s^.  The  target  of  these 
simulations  were  Milky  Way  sized  halos  with  higher  mass 
resolution  than  our  simulation  but  for  a  CDM  cosmology 
only.  There  is  a  large  variation  in  magnitude  between  the 
simulations  that  cannot  be  explained  as  a  resolution  arti¬ 
fact  or  by  differences  in  the  adopted  cosmological  param¬ 
eters.  Diemand  et  al.  argued  the  difference  between 
Via  Lactea  I  and  H  lies  within  the  halo-to-halo  varance. 
The  Aquarius  simulation  suite  included  6  unique  halos 
with  a  rms  scatter  of  only  10%  about  the  mean.  The 
abundance  of  satellites  in  Via  Lactea  I  is  69%  lower  than 
the  Aquarius  mean  and  Via  Lactea  H  is  31%  lower  which 
led  Springel  et  al.  [s^  to  argue  that  the  difference  be¬ 
tween  Aquarius  and  Via  Lactea  is  due  to  a  systematic 
error  in  the  Via  Lactea  realization  of  the  initial  condi¬ 
tions.  Our  CDM  simulations  have  satellite  abundances 
between  Via  Lactea  I  &  H.  Both  Via  Lactea  and  our  sim¬ 
ulations  used  initial  conditions  generated  with  GRAFIC2 
which  may  be  why  our  results  disagree  with  Aquarius. 

On  the  other  hand,  Ishiyama  et  al.  examined  the 
satellite  abundances  of  68  simulated  Milky  Way  sized 
halos.  They  claim  their  simulations  agree  with  both 
Via  Lactea  and  Aquarius.  They  also  found  a  correla¬ 
tion  between  satellite  abundance  and  halo  concentration 
and  between  the  maximum  value  of  the  radius  enclosing 
half  of  the  final  virial  mass.  This  suggests  the  number 
of  satellites  is  determined  by  the  formation  epoch  of  the 
halo  and  at  least  partly  by  the  assembly  history  of  the 
halo.  It  is  possible  the  small  scatter  in  the  Aquarius 
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TABLE  I.  Properties  of  simulated  Milky  Way  halos. 


Simulation 

M200 

R200 

M50 

R50 

aso 

'^max 

N200 

N50 

[kpc] 

[lO^^Mo] 

[kpc] 

[km/s]  [km/s] 

Set  A 

CDM  lo 

1.2919 

218.49 

1.6792 

378.49 

138.14 

183.34 

1,  760,  904 

2,  288,  747 

5  keV  lo 

1.2991 

218.90 

1.6829 

378.77 

138.24 

184.01 

1,  770,  649 

2,  293,  847 

4  keV  lo 

1.2937 

218.63 

1.6838 

378.90 

138.25 

184.62 

1,  763,  346 

2,295,103 

3  keV  lo 

1.2882 

218.36 

1.6855 

379.04 

138.29 

184.21 

1,755,895 

2,297,320 

2  keV  lo 

1.3115 

219.59 

1.6703 

377.81 

137.89 

182.55 

1,  787,  535 

2,  276,  577 

1  keV  lo 

1.3056 

219.32 

1.6616 

377.26 

137.64 

180.94 

1,779,518 

2,  264,  844 

CDM  hi 

1.5278 

231.10 

2.0354 

403.56 

147.28 

192.27 

16,  659,  568 

22, 194, 175 

4  keV  hi 

1.5168 

230.55 

2.0403 

403.97 

147.39 

190.74 

16,  538,  776 

22,247,815 

2  keV  hi 

1.5272 

231.10 

2.0299 

403.29 

147.14 

184.35 

16,652,697  22,134,513 

1  keV  hi 

1.5493 

232.19 

2.0254 

402.88 

147.05 

179.96 

16,  893,  692 

22,  084,  649 

Set  B 


CDM  lo 

1.6440 

236.85 

2.1324 

409.86 

149.59 

196.54 

2,  240,  809 

2,  906,  506 

5  keV  lo 

1.6424 

236.71 

2.1255 

409.45 

149.42 

196.49 

2,  238,  592 

2,897,056 

4  keV  lo 

1.6417 

236.71 

2.1213 

409.18 

149.33 

196.82 

2,  237,  605 

2,891,353 

3  keV  lo 

1.6358 

236.44 

2.1186 

409.04 

149.25 

196.40 

2,  229,  660 

2,  887,  665 

2  keV  lo 

1.6128 

235.34 

2.0935 

407.40 

148.67 

196.15 

2, 198,  289 

2,  853,  496 

1  keV  lo 

1.5798 

233.70 

2.0750 

406.16 

148.23 

193.26 

2, 153,  349 

2,828,261 

CDM  hi 

1.5155 

230.41 

1.9966 

400.96 

146.35 

194.99 

16,  525, 148 

21,771,052 

4  keV  hi 

1.4883 

229.04 

1.9819 

400.00 

145.98 

189.56 

16,228,  561 

21,611,248 

2  keV  hi 

1.4759 

228.49 

1.9667 

399.04 

145.59 

186.01 

16,  093,  504 

21,444,  730 

1  keV  hi 

1.3843 

223.56 

1.8656 

392.05 

143.06 

180.09 

15,094,  535 

20,  342, 179 

halos  may  be  a  selection  effect.  The  6  Aquarius  halos 
were  chosen  because  they  hosted  normal  spiral  galaxies  in 
semi-analytic  modeling  and  had  no  massive  neighbors  at 
z  =  0,  suggesting  the  halos  had  similar  formation  epochs 
and  accretion  histories.  Ishiyama’s  halos  were  not  se¬ 
lected  by  environment  or  accretion  history.  Indeed,  their 
halos  with  the  highest  abundance  of  satellites  have  late 
formation  epochs  and  are  still  undergoing  major  merg¬ 
ers,  but  such  halos  are  unlikely  to  be  representative  of 
the  Milky  Way.  Satellite  mergers  heat  and  potentially 
destroy  the  disks  of  spiral  galaxies  [H,  [s^.  The  thick 
disk  of  the  Milky  Way  is  thought  to  have  been  created  by 
a  merger  event  about  12  billion  years  ago.  The  10  —  12 
billion  year  old  thin  disk  indicates  the  Milky  Way  has 
not  suffered  a  major  merger  since  2^1.  Aquarius’  semi- 
analytic  modeling  shows  their  halos  host  Milky  Way-type 
spiral  galaxies  and  should  be  characteristic  of  the  Milky 
Way.  The  Aquarius  simulations  also  have  the  highest 
abundances  of  simulations  attempting  to  represent  the 
Milky  Way.  If  the  Aquarius  results  are  correct  and  their 
satellite  abundances  represent  the  Milky  Way,  our  sim¬ 
ulations  can  be  brought  into  agreement  by  multiplying 
the  set  A  abundances  by  2.5  and  set  B  by  2.0. 


B.  Poisson  Halos 

We  discussed  in  Section  U  that  the  discrete  sampling  of 
the  matter  density  field  in  N-body  simulations  introduces 
Poisson  noise  to  the  power  spectrum  of  density  perturba¬ 
tions.  The  Poisson  noise  results  in  the  formation  of  arti¬ 
ficial  small  mass  halos  in  WDM  cosmologies.  Plotted  in 
Figure  [7]  are  the  circular  velocity  functions  compared  to 
analytical  predictions  using  Press-Schechter  theory  (40l|. 
Barkana,  Haiman,  &  Ostriker  [4l|  showed  the  effects  of 
WDM  on  Press  Schechter  theory.  We  extend  their  work 
by  also  including  the  Poisson  noise  from  the  discreteness 
of  the  simulations.  Our  model  of  the  effect  of  Poisson 
noise  on  the  power  spectrum  is  analogous  to  the  one  in 
Afshordi,  McDonald,  &  Spergel  [i^ ,  that  considered  the 
case  of  dark  matter  composed  of  primordial  black  holes. 
They  showed  that  the  fluctuation  in  the  number  of  pri¬ 
mordial  black  holes  due  to  Poisson  noise  leads  to  the 
addition  of  a  constant  term  to  the  linear  power  spec¬ 
trum  that  becomes  dominant  at  small  scales.  We  add 
a  constant  Poisson  noise  term  to  the  power  spectrum  in 
our  Press-Schechter  calculations.  The  addition  of  Poisson 
noise  produces  the  upturns  at  low  velocities  in  Figure  [T] 
These  upturns  are  also  apparent  in  the  simulations  of 
Bode,  Ostriker,  &  Turok  and  Barkana,  Haiman,  & 
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FIG.  4.  Portraits  of  the  simulated  Milky  Way  halo  at  2:  =  0  in  the  set  B  low  resolution  simulations.  Images  are  1  Mpc  on  a 
side. 


Ostriker  [4l|  but  were  attributed  to  the  fragmentation  of 
matter  filaments.  The  upturns  occur  at  velocities  that 
depend  on  the  cosmology.  For  our  warmest  cosmology 
Poisson  halos  are  important  only  at  circular  velocities 
less  than  about  6  km/s. 

The  vertical  lines  in  Figure  0  show  where  Vmax  = 


8  km/s.  Below  this  scale  the  high  resolution  CDM 
simulations  begin  to  fall  away  from  the  Press-Schecter 
lines  due  to  the  resolution  limits  of  the  simulation.  For 
Vmax  >  8  km/s  our  simulations  are  reasonably  complete 
within  i?5o  of  each  Milky  Way  although  numerical  de¬ 
struction  of  a  small  fraction  of  satellites  in  the  inner 
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FIG.  5.  Milky  Way  halo  and  satellites  in  the  set  B  high  reso¬ 
lution  CDM  simulation  before  (top)  and  after  (bottom)  sub¬ 
tracting  the  Milky  Way  profile.  Only  the  satellites  remain 
after  subtraction.  The  halo  was  subtracted  to  a  distance  of 
500  kpc  from  the  Milky  Way  center. 


Milky  Way  wouldn’t  be  apparent  in  Figure  [71  especially 
for  the  CDM  and  4  keV  cosmologies.  Before  compar¬ 
ing  our  simulations  to  observations  of  the  Milky  Way  we 
need  to  determine  to  what  distance  our  simulations  are 
convergent. 


C.  Convergence  Study 

Satellites  orbiting  in  the  halo  of  a  larger  galaxy  are  de¬ 
stroyed  by  tidal  stripping  and  heating  through  encounters 
with  other  satellites.  Satellites  in  simulations  are  also  de¬ 
stroyed  artificially  by  numerical  effects  that  become  dom¬ 
inant  for  poorly  resolved  halos  in  the  inner  halo  region. 
There  will  therefore  be  a  radius  inside  of  which  our  sim¬ 
ulations  will  not  converge  to  a  realistic  representation  of 
the  actual  Milky  Way. 

To  determine  the  convergence  of  our  simulations  and 
have  an  idea  of  the  variance  of  the  results,  we  performed 
two  simulations  at  lower  and  higher  resolution  and  we 
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FIG.  6.  Gumulative  mass  functions  for  satellites  in  our  low 
(dotted)  and  high  (solid)  resolution  set  B  (top)  and  set  A  (bot¬ 
tom)  simulations.  Satellite  masses  are  characterized  by  their 
circular  velocities  and  have  been  normalized  to  the  circular 
velocity  of  the  host  halo  at  a  radius  enclosing  an  overden¬ 
sity  of  50.  Straight  sloped  lines  are  the  fitting  formula  from 
Via  Lactea  I,  Via  Lactea  II,  and  the  Aquarius  simulations 
(left  to  right). 


also  simulated  two  different  realizations  of  a  Milky  Way¬ 
sized  galaxy.  We  performed  convergence  studies  follow¬ 
ing  the  argument  elucidated  below,  in  combination  with 
results  of  published  high  resolution  simulations  found  in 
the  literature.  Using  the  work  of  Moore  et  ah,  De  Lu¬ 
cia  et  ah,  and  Ishiyama  et  al.  d,  113,  Ell  we  will  assume 
that  the  shape  of  the  cumulative  satellite  mass  function 
for  host  halos  of  different  masses  is  nearly  constant  and 
the  total  number  of  satellites  scales  linearly  with  the  host 
mass.  If  the  simulations  are  convergent  the  cumulative 
circular  velocity  function  for  satellites,  N(R),  within  a 
given  Galactocentric  radius,  R,  should  be  proportional 
to  the  enclosed  mass,  M(R),  and  a  function  of  R  that 
represents  the  fraction  of  satellites  that  survive  destruc¬ 
tion  from  physical  effects: 

N(R)  =  f(R)M(R),  (9) 

where  f(R)  oc  i?“.  The  normalization  of  f(R)  can  be  set 
using  values  of  N(R)  and  M(R)  at  a  distance  Rq: 

(j^) 

The  mass  functions  normalized  in  this  way  will  be  con¬ 
stant  with  radius  where  the  simulations  are  convergent. 
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FIG.  7.  Cumulative  mass  functions  for  satellites  in  our  low 
(dotted)  and  high  (solid)  resolution  set  B  (top)  and  set  A 
(bottom)  simulations.  Mass  functions  calculated  from  Press- 
Schechter  theory  shifted  vertically  to  match  the  CDM  curves 
have  been  added.  The  WDM  simulations  are  well  fit  by  the 
theory.  The  addition  of  a  constant  power  Poisson  noise  creates 
the  upturns  in  the  WDM  cosmologies  at  velocites  <  6  km/s. 
The  level  of  Poisson  noise  was  adjusted  to  match  the  upturn 
at  1  keV  in  set  B  but  also  matches  the  2  keV  well  and  the  set 
A  simulations.  Below  8  km/s  (vertical  line)  the  CDM  lines  fall 
below  the  Press-Schechter  calculations  due  to  the  resolution 
limits  of  the  simulations. 


Where  numerical  effects  destroy  satellites  the  mass  func¬ 
tions  will  normalize  to  a  lower  value.  We  expect  that 
a  is  constant  because  there  isn’t  any  characteristic  scale 
for  the  destruction  rate  in  dark  matter  only  simulations. 
Hence,  a  can  be  determined  at  large  radii  where  conver¬ 
gence  is  certain. 

Figure  [8]  shows  the  normalized  mass  functions  for  the 
simulations.  The  normalization  constants  Mo,Rq  have 
been  chosen  at  200  kpc  and  the  value  of  a  (0.55)  was 
adjusted  by  hand  until  a  good  fit  was  achieved  for  the 
set  B  mass  functions  above  200  kpc  in  the  high  reso¬ 
lution  CDM  cosmology  at  circular  velocities  >  8  km/s 
(vertical  lines).  This  a  also  provides  a  good  fit  for  the 
WDM  cosmologies  and  for  the  set  A  simulations,  al¬ 
though  the  1  and  2  keV  mass  functions  have  a  wider 
scatter  due  to  the  smaller  numbers  of  satellites  in  these 
simulations.  The  mwDM  =  4  keV  simulation  is  conver¬ 
gent  for  Vmax  >  8  km/s  to  distances  >  100  kpc.  At 
75  kpc  the  effects  of  numerical  resolution  are  apparent. 
The  same  value  of  a  has  been  used  in  the  normalization 
of  the  low  resolution  sets  and  appears  to  provide  a  good 


fit  for  the  mass  functions  >  250  kpc  (thin  solid  lines). 
The  effects  of  numerical  resolution  on  the  destruction  of 
satellites  are  apparent  at  larger  distances  in  these  simu¬ 
lations:  <  200  kpc  for  CDM  and  <  150  kpc  for  WDM. 


D.  Comparison  to  Observations 

Before  the  Sloan  Digital  Sky  Survey  there  were  only 
12  classically  known  satellite  galaxies  to  the  Milky  Way. 
Sixteen  new  satellites  have  been  discovered  in  the  SDSS, 
currently  in  Data  Release  7.  We  list  all  known  Milky 
Way  satellites  in  Table  HIl  We  use  the  satellite  distances 
given  in  Table UTl as  their  Galactocentric  distances.  When 
comparing  the  observed  satellites  to  the  simulations  we 
must  correct  the  SDSS  dwarfs  for  completeness.  The 
primary  incompleteness  of  the  SDSS  is  its  sky  coverage, 
28.3%  (11663  deg^).  Second,  being  a  magnitude  limited 
survey,  the  SDSS  has  a  luminosity  bias.  The  detection  ef¬ 
ficiency  of  dwarfs  in  the  SDSS  is  a  function  of  dwarf  size, 
luminosity,  distance,  and  Galactic  latitude  as  shown  by 
Walsh  et  al.  HL  An  approximate  expression  is  given  in 
Tollerud  et  al.  [45j|  (using  the  work  of  Koposov  et  al.  0) 
for  the  distance  which  galaxies  of  luminosity  >  L  are 
completely  detected:  d  «  66kpc(L/1000  Galax¬ 

ies  with  L  >  10^  Lq  should  be  approximately  complete 
to  200  kpc,  with  L  >  2300  Lq  to  100  kpc.  The  distance 
range  100—200  kpc  is  thus  suited  for  comparisons  because 
our  simulations  are  convergent  and  the  observations  are 
nearly,  but  not  quite,  complete.  For  our  analysis  we  will 
only  use  satellites  with  distances  <  200  kpc.  For  compar¬ 
ison  with  the  simulations  we  do  a  first  order  correction 
for  the  sky  coverage  of  the  SDSS  assuming  an  isotropic 
distribution  of  satellites.  The  SDSS  covers  0.283  of  the 
sky  so  for  every  satellite  detected  we  assume  there  are  ac¬ 
tually  3.54  satellites  at  that  distance.  Gombined  with  the 
classic  Milky  Way  satellites  this  forms  our  observed  data 
set.  We  also  implemented  a  conservative  luminosity  cor¬ 
rection  for  the  SDSS  dwarfs  using  the  formulas  in  Walsh 
et  al.  [i^ .  This  adds  2  satellites  making  a  total  of  61  ±  13 
satellites  within  200  kpc  where  the  error  only  includes 
Poisson  statistics  of  the  SDSS  dwarfs:  3.54i/Amss.  The 
extra  2  satellites  are  at  distances  150  —  200  kpc  and  do 
not  affect  our  conclusions.  We  note  that  we  get  the  same 
result  if  we  only  include  dwarfs  detected  in  the  Data  Re¬ 
lease  5  footprint  and  correct  for  the  smaller  sky  cover- 
age  (8000  deg^).  We  also  note  the  formulas  in  Walsh  et 
al.  j44|  assume  the  size-luminosity  distribution  of  known 
dwarfs  is  representative  of  all  satellites.  There  may  be  a 
population  of  dwarfs  with  surface  brightnesses  below  the 
detection  limit  of  the  SDSS  [62^0 . 

Willman  1  is  an  exceptional  case  in  that  it  may  not  be  a 
dark  matter  dominated  dwarf  galaxy  but  a  globular  clus¬ 
ter  undergoing  tidal  disruption.  Its  size  and  luminosity 
are  intermediate  between  Milky  Way  dwarfs  and  globular 
clusters  [s^.  Although  it  has  a  large  metallicity  spread 
unlike  the  stellar  population  of  a  globular  cluster,  when 
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FIG.  8.  (Left).  Normalized  mass  functions  for  set  B  high  {top)  and  low  (bottom)  resolution  simulations.  Thin  solid  lines  are 
distances  of  400,  300,  and  250  kpc.  Thick  solid  line  is  200  kpc,  dotted  line  is  150  kpc,  dashed  line  is  100  kpc,  dot-dashed  line  is 
75  kpc.  The  WDM  cosmologies  have  been  shifted  down  vertically  for  clarity.  The  value  a  —  0.55  was  set  by  the  high  resolution 
simulation  and  provides  good  normalization  for  the  low  resolution  as  well  but  the  effects  of  incompleteness  become  apparent 
at  much  larger  radii  (150  —  200  kpc  compared  to  75  —  100  kpc  for  high  res).  (Right).  Same  as  the  left  panel,  but  for  the  set  A 
high  (top)  and  low  (bottom)  resolution  simulations.  The  value  a  =  0.55  also  provides  good  normalization  for  this  set  of  halos. 


deriving  constraints  on  the  dark  matter  particle  mass  we 
will  consider  both  including  and  excluding  Willman  1  as 
a  Milky  Way  satellite. 

When  comparing  observations  and  simulations  we  need 
to  know  the  maximum  circular  velocities  of  the  dark 
matter  halos  the  observed  galaxies  are  presumably  em¬ 
bedded  in.  Ricotti  &  Gnedin  found  in  simulations 
that  the  maximum  circular  velocities  of  satellites  are  at 
least  twice  the  velocity  dispersion  of  the  stellar  compo¬ 
nent.  Assuming  the  stellar  velocity  dispersions  of  the 
observed  dwarfs  are  times  the  line-of-sight  velocity 
dispersions  (<Jstar  in  Table  ITTll .  then  all  dwarfs  with  mea¬ 
sured  velocity  dispersions  have  max  circular  velocities 
greater  than  8  km/s.  We  assume  dwarfs  without  mea¬ 
sured  velocity  dispersions  are  similar  to  the  other  known 
dwarfs  and  conservatively  conclude  all  dwarfs  reside  in 
dark  matter  halos  with  max  circular  velocities  greater 
than  8  km/s.  An  alternative  approach  is  the  work  of 
Wolf  et  al.  relating  the  circular  velocity  at  half  light 
radius  to  a  star-  Vc(r  1/2)  =  VScTstar-  All  the  observed 
dwarfs  except  Leo  V  have  circular  velocities  at  half  light 
radius  about  6  km/s  or  greater.  Since  the  max  circular 
velocity  must  be  greater  than  or  equal  to  the  half  light 
circular  velocity  we  will  also  consider  that  all  observed 
dwarfs  reside  in  halos  with  Vmax  >  6  km/s. 

In  the  left  panel  of  Figure  [9]  we  plot  histograms  of  the 
number  of  satellites  with  distance  for  the  observed  and 


simulation  data  sets  with  an  8  km/s  max  circular  veloc¬ 
ity  cut.  The  upward  arrows  on  the  observed  data  bars 
indicate  these  are  only  lower  limits  due  to  the  surface 
brightness  limits  of  the  SDSS,  it  is  possible  there  are  more 
dwarfs  yet  to  be  discovered.  The  8  km/s  cut  to  the  simu¬ 
lation  data  avoids  Poisson  halos  and  assures  the  high  res¬ 
olution  simulations  are  complete  to  at  least  r  =  100  kpc. 
The  low  resolution  simulations  are  also  plotted  in  these 
figures  but  they  are  complete  only  to  150  kpc.  Focusing 
on  the  100  —  200  kpc  bins  it  is  clear  the  1  keV  has  far 
too  few  satellites  to  match  the  observations.  The  2  keV 
simulations  can  be  consistent  with  the  observations  if  the 
simulations  are  incomplete  below  100  kpc  or  the  sky  cor¬ 
rection  has  overestimated  the  number  of  satellites  in  the 
inner  Milky  Way.  The  4  keV  simulations  can  also  be 
consistent  with  the  observations  although  they  may  re¬ 
quire  some  of  the  dark  matter  halos  to  not  host  luminous 
galaxies.  It  is  not  clear  from  this  plot  how  variance  in 
the  satellite  abundances  for  the  simulated  halos  may  af¬ 
fect  the  results. 

The  number  of  satellites  in  the  simulations  can  be  cor¬ 
rected  for  completeness  using  the  convergence  equation, 
Eq.  (fTUl).  We  used  the  mass  and  number  of  satellites 
inside  A50  for  the  normalization  and  calculated  the  num¬ 
ber  of  satellites  in  50  kpc  bins  for  the  high  resolution 
simulations.  The  results  are  shown  in  the  right  panel  of 
Figure  O  The  results  are  very  similar  across  cosmolo- 
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FIG.  9.  (Left)  Number  of  satellites  with  distance  from  the  Milky  Way  grouped  into  50  kpc  bins.  The  simulation  satellites 
have  been  cut  by  circular  velocity  >  8  km/s.  Solid  lines  show  the  data  from  the  high  resolution  and  dashed  lines  show  the  low 
resolution  simulations,  thick  lines  are  the  set  B  simulations.  The  bars  with  arrows  are  the  observed  satellites  after  correcting 
for  the  sky  coverage  of  the  SDSS.  Observations  are  incomplete  at  distances  greater  than  about  50  —  100  kpc  (depending  on 
the  lumiosity  and  surface  brightness  of  the  dwarf),  while  simulations  have  not  converged  for  less  than  about  100  kpc  for  high 
resolution  and  150  kpc  for  low  resolution.  (Right)  Number  of  satellites  with  distance  from  the  Milky  Way  like  the  plot  at  left 
but  calculated  from  the  convergence  equation  lEo.  IIOII.  The  convergence  correction  affects  mainly  the  4  keV  0  —  50  kpc  bin. 


gies,  a  nearly  constant  number  of  satellites  per  bin  from 
50  —  200  kpc  with  about  half  as  many  in  the  0  —  50  kpc 
bin.  The  plots  are  also  very  similar  to  the  simulation 
data  plots  except  for  the  0  —  50  kpc  bin  in  the  4  keV 
cosmology  where  about  10  satellites  were  destroyed  by 
numerical  effects.  The  0  —  50  kpc  bin  is  most  important 
for  constraining  the  dark  matter  particle  mass  because 
the  observations  are  most  complete  in  this  bin. 

In  Figure  (TU]  we  plot  the  number  of  satellites  in  the 
0  —  50  kpc  bin  calculated  from  the  convergence  equa¬ 
tion  for  the  high  resolution  set  B  Milky  Way  as  a  func¬ 
tion  of  the  particle  mass.  To  account  for  the  variance  in 
satellite  abundances  we  also  scale  to  the  Aquarius  and 
Via  Lactea  I  simulations  (shaded  region).  We  assume 
the  Aquarius  abundances  represent  the  upper  limit  for 
the  Milky  Way  since  these  simulations  have  the  largest 
abundances  of  published  Milky  Way-sized  dark  matter 
halos  that  also  host  spiral  galaxies  or  are  isolated  and 
therefore  likely  to  host  spirals.  Some  of  the  halos  in 
Ishiyama  et  al.  [sj  have  higher  abundances  but  they  are 
also  undergoing  major  mergers  and  not  representative 
of  the  Milky  Way.  The  solid  horizontal  line  shows  the 
number  of  observed  satellites  after  applying  sky  coverage 
corrections  to  the  SDSS  data  while  the  dashed  horizon¬ 
tal  lines  show  the  68%  and  95%  confidence  lower  lim¬ 
its  and  the  raw  number  of  observed  satellites  without 
any  correction.  From  the  corrected  observations  we  get 
mwDM  >  3.7  keV  with  >  3.2  keV  at  68%  confidence  and 
>  2.6  keV  at  95%  confidence.  The  raw  observed  line  is  an 
extremely  conservative  limit  as  the  number  of  Milky  Way 
satellites  within  50  kpc  is  almost  certainly  greater  than 
7,  but  it  sets  an  absolute  lower  limit  of  mwDM  >  2  keV. 


FIG.  10.  Number  of  satellites  from  0-50  kpc  calculated  from 
the  convergence  equation  in  the  set  B  high  resolution  simula¬ 
tions  for  Vmax  >  8  km/s  (thick  line  with  points)  compared  to 
the  observed  number  of  Milky  Way  satellites,  including  Will- 
man  1,  with  corrections  (solid  horizontal  line)  and  without 
(“raw”)  and  at  68%  and  95%  confidence  (dashed  horizontal 
lines).  The  gray  shaded  area  is  the  simulation  data  scaled  to 
the  satellite  abundances  of  the  Aquarius  simulations  at  the 
high  end  and  Via  Lactea  I  on  the  low  end.  The  uncorrected 
observations  intersect  the  shaded  region  at  mwDM  ~  2  keV, 
this  is  an  extremely  conservative  lower  limit  on  the  dark  mat¬ 
ter  particle  mass. 

We  repeated  the  same  analysis  using  a  6  km/s  max  cir¬ 
cular  velocity  cut  to  the  simulation  data.  We  also  con¬ 
sidered  the  effects  when  Willman  1  was  excluded  from 
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TABLE  II.  Summary  of  known  Milky  Way  satellites. 


Name  dist  Cstar  Mv  References 

[kpc]  [km/s] 


Classical  (pre-SDSS) 


Sagittar 

24  ±2 

11.4  ±0.7 

-13.4 

[47] 

LMC 

49  ±2 

- 

-18.4 

[47] 

SMC 

58  ±2 

- 

-17.0 

[47] 

Ursa  Minor 

66  ±3 

9.3  ±  1.8 

-8.9 

m 

Draco 

79  ±4 

9.5  ±1.6 

-8.8 

[47] 

Sculptor 

79  ±4 

6.6  ±0.7 

-11.1 

[47] 

Sextans 

86  ±4 

6.6  ±0.7 

-9.5 

m 

Carina 

94  ±5 

6.8  ±  1.6 

-9.3 

[47] 

Fornax 

138  ±8 

10.5  ±  1.5 

-13.2 

[47] 

Leo  II 

205  ±  12 

6.7±  1.1 

-9.6 

m 

Leo  I 

270  ±  30 

8.8  ±0.9 

-11.9 

[47] 

Phoenix 

405  ±  15 

- 

-10.1 

[47] 

SDSS  discovered 

Segue  I 

23  ±2 

4.3  ±  1.2 

-1.5 

[48] 

Ursa  Major  II 

30  ±5 

6.7  ±  1.4 

-3.8 

[49,^ 

Segue  II 

~  35 

3.4  ±2.0 

-2.5 

[^ 

Willman  I 

38  ±7 

4.3/);3 

-2.5 

[49,^ 

Coma  Berenics 

44  ±4 

4.6  ±0.8 

-3.7 

[OT,^ 

Bootes  II 

60±  10 

- 

-3.1 

Mi 

Bootes  I 

62  ±3 

6-5/?:4 

-5.8 

[49] 

Pisces  I 

80±  14 

- 

- 

MM 

Ursa  Major  I 

106/® 

7.6  ±  1.0 

-5.6 

[^ 

Hercules 

140/^ 

5.1  ±0.9 

-6.0 

[OT,^ 

Canes  Venatici  II 

150/)^ 

4.6  ±  1.0 

-4.8 

[OT,^ 

Leo  IV 

160/)® 

3.3  ±  1.7 

-5.8 

[50,^ 

Leo  V 

175  ±9 

2.4  ±  1.8 

-5.2 

MM 

Pisces  II 

~  180 

- 

-5.0 

[^ 

Canes  Venatici  I 

220/)^ 

7.6  ±0.4 

-7.9 

[M,  60] 

Leo  T 

~  420 

7.5  ±1.6 

-7.1 

the  observed  data  set  for  both  the  6  km/s  and  8  km/s 
analysis.  We  present  the  results  in  Table  IIIII  In  the 


TABLE  III.  Dark  matter  particle  mass  constraints  (in  keV) 
for  simulated  halo  max  circular  velocity  cuts  of  8  and  6  km/s 
and  including  or  excluding  Willman  1  from  the  observed  data 
set. 


Willman  1? 

Vmax  >  8  km/s 
Included  Excluded 

Vmax  >  6  km/s 

Included  Excluded 

raw  observed 

>  2.0 

>  1.9 

>  1.8 

>  1.6 

95%  confidence 

>  2.6 

>  2.3 

>  2.3 

>  2.1 

68%  confidence 

>  3.2 

>  2.6 

>  2.6 

>  2.3 

sky  corrected 

>  3.7 

>  3.2 

>  3.0 

>  2.6 

most  conservative  case,  where  Willman  1  is  not  a  dark 
matter  dominated  dwarf  galaxy  and  all  observed  satel¬ 
lites  correspond  to  dark  matter  halos  with  max  circular 
velocities  >  6  km/s,  we  can  say  mwDM  >  2.1  keV  with 
95%  confidence.  We  adopt  this  as  our  formal  limit  for 
this  work. 


IV.  DISCUSSION 

We  found  in  the  preceding  section  a  conservative  lower 
limit  of  my/ DM  >  2.1  keV  on  the  dark  matter  particle 
mass.  We  also  found  the  1  keV  WDM  simulation  has  too 
few  satellites  to  match  the  Milky  Way  observations.  This 
agrees  with  the  semi-analytic  modeling  and  Milky  Way 
satellite  luminosity  functions  in  WDM  cosmologies  work 
of  Maccio  &  Fontanot  0 . 

Our  result  can  also  be  compared  to  limits  on  the 
particle  mass  from  the  Lyman-a  forest  in  high  redshift 
quasars.  Lyman-a  absorption  by  neutral  hydrogen  along 
the  line  of  sight  to  distant  quasars  over  redshifts  2-6 
probes  the  matter  power  spectrum  in  the  mildly  nonlin¬ 
ear  regime  on  scales  1-80  Mpc/h.  Viel  et  al.  [lo,  [63, 
have  numerically  modeled  the  Lyman-a  forest  flux  power 
spectra  for  varied  cosmological  parameters  and  compared 
to  observed  quasar  forests  to  obtain  lower  limits  on  the 
dark  matter  particle  mass.  Their  2006  work  used  low 
resolution  spectra  for  3035  quasars  (2.2  <  z  <  4.2)  from 
the  SDSS  [6^  and  found  a  2cr  lower  limit  of  2  keV  for 
a  thermal  WDM  particle.  This  limit  agrees  with  our 
results  that  a  2.1  keV  particle  is  the  lower  limit  that 
can  reproduce  the  observed  number  of  Milky  Way  satel¬ 
lites.  However,  their  latest  work  uses  high  resolution 
spectra  for  55  quasars  (2.0  <  z  <  6.4)  from  the  Keck 
HIRES  spectrograph  in  addition  to  the  SDSS  quasars. 
With  the  new  data  they  report  a  lower  limit  of  4  keV 
{2a).  A  caveat  arises  in  Viel  et  al.  [t^  who  show  the  flux 
power  spectrum  from  the  SDSS  data  prefer  larger  values 
of  the  intergalactic  medium  (IGM)  temperature  at  mean 
density  than  expected  from  photoionization.  The  flux 
power  spectrum  temperature  is  also  higher  than  that  de¬ 
rived  from  an  analysis  of  the  flux  probability  distribution 
function  of  18  high  resolution  spectra  from  UVES/VLT 
and  also  higher  than  constraints  from  the  widths  of  ther¬ 
mally  broadened  absorption  lines  [zllzl-  This  could  be 
explained  by  an  unaccounted  for  systematic  error  in  the 
SDSS  flux  power  spectrum  data  which  may  also  affect 
the  derived  dark  matter  particle  mass  limits. 

Using  the  scaling  relation  for  sterile  neutrinos  we  find 
a  lower  limit  TOs  >  11.8  keV  with  95%  confidence  for  a 
DW  produced  sterile  neutrino  particle.  Scaling  to  the 
other  production  mechanisms  we  get  m-s  >  7.9  keV  for 
the  SE  mechanism  and  rUg  >  2.6  keV  for  Higgs  decay 
sterile  neutrinos.  Sterile  neutrinos  are  expected  to  ra- 
diatively  decay  to  a  photon  and  a  lighter  mass  neutrino. 
This  decay  should  produce  an  X-ray  spectral  line  with 
energy  =  ms/2.  Abazajian  et  al.  |73i|  used  Chandra 
spectra  of  the  the  unresolved  component  of  the  cosmic  X- 
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ray  background  in  the  Chandra  Deep  Fields  North  and 
South  to  search  for  a  contribution  from  the  decay  of  ster¬ 
ile  neutrinos  in  the  dark  matter  halo  of  the  Milky  Way. 
Their  non-detection  places  an  upper  limit  on  the  sterile 
neutrino  mass  that  depends  on  the  model  for  the  Milky 
Way  halo.  Using  their  lowest  mass  models  they  adopt 
a  conservative  limit  of  mg  <  5.7  keV  (2a).  Our  lower 
limit  for  DW  and  SF  produced  neutrinos  is  significantly 
above  the  X-ray  upper  limit  implying  the  dark  matter  is 
not  made  entirely  of  sterile  neutrinos  produced  by  these 
mechanisms.  If  sterile  neutrinos  compose  only  part  of  the 
dark  matter  the  X-ray  upper  mass  limit  can  be  higher  but 
our  lower  mass  limits  do  not  necessarily  apply  to  mixed 
dark  matter  cosmologies.  Higgs  decay  produced  sterile 
neutrinos  can,  however,  constitute  all  the  dark  matter 
and  be  consistent  with  the  X-ray  and  Milky  Way  satel¬ 
lite  constraints.  These  conclusions  agree  with  the  results 
of  Suzaku  observations  of  Ursa  Minor  [t^  . 


V.  SUMMARY 

We  have  conducted  N-body  simulations  of  the  forma¬ 
tion  of  a  Milky  Way-sized  dark  matter  halo  in  CDM  and 
WDM  cosmologies.  Such  simulations  are  complicated  by 
the  formation  of  nonphysical  small  mass  halos  due  to  the 
discreteness  of  the  initial  conditions.  We  have  shown  that 
these  halos  can  be  modeled  by  including  Poisson  noise  in 
Press-Schechter  calculations  and  some  do  survive  to  end 
up  inside  the  halo  of  the  Milky  Way  but  with  sufficient 
resolution  they  are  only  important  at  small  scales  and 
can  be  avoided  with  an  appropriate  circular  velocity  or 
mass  cut. 

We  have  studied  the  number  of  satellite  halos  as  a  func¬ 
tion  of  distance  from  the  Milky  Way.  The  4  keV  and 
2  keV  WDM  simulations  can  adequately  reproduce  the 
observed  number  of  satellites  at  hundreds  of  kiloparsecs 
while  the  1  keV  is  severely  dehcient.  Correcting  our  sim¬ 
ulations  for  numerical  destruction  in  the  inner  50  kpc  and 
scaling  to  the  higher  satellite  abundances  of  the  Aquar¬ 
ius  simulations,  we  find  a  very  conservative  lower  limit  of 
>  2.1  keV  at  95%  confidence.  This  agrees  with  the  ear¬ 
lier  Lyman-a  forest  modeling  work  of  Viel  et  al.  that 
mwDM  >  2  keV  but  the  two  methods  are  independent 
and  almost  certainly  are  subject  to  different  systematic 
errors  if  any  exist.  Their  latest  work  raises  the  limit 
to  mwDM  >  4  keV  but  problems  with  the  derived  IGM 
temperature  and  mean  density  may  indicate  the  SDSS 
spectra  they  used  suffer  a  systematic  error  0. 


Our  lower  limit  of  2.1  keV  for  a  thermal  dark  matter 
particle  scales  to  lower  limits  of  11.8,  7.9,  2.6  keV  for  DW, 
SF,  and  Higgs  decay  produced  sterile  neutrinos.  Sterile 
neutrinos,  if  they  exist,  are  expected  to  decay  into  X- 
rays  and  active  neutrinos.  Observations  of  the  unresolved 
cosmic  X-ray  background  combined  with  models  of  the 
Milky  Way  halo  set  an  upper  limit  of  5.7  keV  on  the 
mass  of  the  sterile  neutrino  |73l| ,  assuming  they  constitute 
100%  of  the  dark  matter.  Sterile  neutrinos  are  viable 
dark  matter  candidates  but  if  they  are  produced  by  the 
DW  or  SF  mechanisms  they  cannot  account  for  all  the 
dark  matter. 

Our  simulations  followed  the  formation  of  two  Milky 
Way-sized  halos  in  CDM  and  WDM  cosmologies.  Nu¬ 
merical  simulations  show  significant  variance  (100%)  in 
the  number  of  satellites  in  Milky  Way-sized  halos,  an  ef¬ 
fect  that  can  be  easily  quantified  using  published  studies 
and  incorporated  in  our  results.  To  account  for  the  effect 
of  variance  we  counted  satellites  in  50  kpc  bins  and  scaled 
our  abundances  to  the  larger  values  of  the  Aquarius  sim¬ 
ulation  suite  which  were  selected,  using  semi-analytic 
modeling,  to  host  Milky  Way-type  spiral  galaxies  and 
showed  <  10%  variance  in  their  number  of  satellites.  Our 
constraint  is  a  conservative  lower  limit  since  we  only  cor¬ 
rect  the  number  of  SDSS  dwarfs  for  sky  completeness.  An 
analysis  that  takes  into  account  the  surface  brightness 
limits  of  the  observational  data  may  allow  tighter  con¬ 
straints  however  the  analysis  would  be  somewhat  model 
dependent. 

We  have  demonstrated  how  N-body  simulations  of 
the  Milky  Way  and  its  satellites  can  set  limits  on  the 
dark  matter  particle  mass  comparable  to  complementary 
methods  such  as  modeling  the  Lyman-a  forest.  These 
limits  are  helped  greatly  by  the  discovery  of  many  new 
Milky  Way  satellites  in  the  SDSS.  There  may  still  be 
a  population  of  low  luminosity,  low  surface  brightness 
dwarf  galaxies  undetectable  by  the  SDSS  Future 

surveys  like  LSST,  DES,  PanSTARRS,  and  SkyMapper 
have  the  potential  to  discover  many  more  Milky  Way 
satellites  and  further  improve  constraints  on  the  mass  of 
the  dark  matter  particle. 
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Appendix 

We  used  the  BBKS  formula  for  the  CDM  transfer 
function  when  generating  the  initial  conditions  for  our 
simulations.  This  formula  assumes  a  baryon  density  of 
zero.  Eisenstein  &  Hu  [l^  calculated  transfer  functions 
for  CDM  cosmologies  that  include  baryon  physics.  We 
plot  the  power  spectra  for  the  fitting  formula  of  Eisen¬ 
stein  &  Hu  (with  fib  =  0.04)  and  BBKS  in  the  left  panel 
of  Figure  [TT]  which  shows  that,  for  a  fixed  value  of  erg, 
BBKS  underestimates  power  on  scales  fc  ^  0.1  but  the 
power  spectra  are  nearly  identical  for  scales  ^  14  Mpc. 
With  flm  =  0.238  a  Milky  Way-sized  halo  with  mass 
^  2  X  10^^  Mq  would  form  from  a  spherical  region  with 
diameter  4.8  Mpc  (fc  =  0.28  h/Mpc).  This  scale  is  well 
in  the  region  where  the  power  spectra  are  nearly  equal 
and  we  would  not  expect  the  use  of  BBKS  to  bias  our 
results. 

To  check  if  BBKS  might  affect  the  number  of  satellites 
we  reran  the  set  B  low  resolution  CDM,  1  keV  WDM, 
2  keV  WDM,  and  high  resolution  CDM  simulations  us¬ 
ing  initial  conditions  generated  from  the  formula  of  Eisen¬ 
stein  &  Hu.  The  right  panel  of  Figure  [TT]  compares  the 
mass  function  of  satellites.  The  results  are  in  good  agree¬ 
ment  and  we  conclude  that  the  use  of  BBKS  has  not 
introduced  a  systematic  error  into  our  results. 
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FIG.  11.  Left  Comparison  of  CDM  power  spectra  calculated  from  the  fitting  formula  of  BBKS  and  Eisenstein  &  Hu  (EH97). 
On  scales  k  >  0.1  h/Mpc  (<  14  Mpc)  the  power  spectra  are  nearly  identical.  The  vertical  line  is  the  diameter  of  a  spherical 
region  with  density  Qmpc  enclosing  a  Milky  Way-sized  mass  2  x  10^^  Mq  5  Mpc).  This  scale  is  well  within  the  range  where 
the  power  spectra  are  nearly  equal.  (Right)  Mass  function  comparison  for  set  B  simulations  using  fitting  formula  for  CDM 
transfer  function  from  BBKS  (dotted)  and  Eisenstein  &  Hu  (solid).  The  results  are  in  good  agreement 


